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Abstract. We consider collisions of particles advected in a fluid. As already pointed 
out by Smoluchowski [Z. f. physik. Chemie, XCII, 129-168, (1917)], macroscopic 
motion of the fluid can significantly enhance the frequency of collisions between the 
suspended particles. This effect was invoked by Saffman and Turner [J. Fluid Mech., 
1, 16-30, (1956)] to estimate collision rates of small water droplets in turbulent rain 
clouds, the macroscopic motion being caused by turbulence. Here we show that the 
Saffman- Turner theory is unsatisfactory because it describes an initial transient only. 
The reason for this failure is that the local flow in the vicinity of a particle is treated 
as if it were a steady hyperbolic flow, whereas in reality it must fluctuate. We derive 
exact expressions for the steady-state collision rate for particles suspended in rapidly 
fluctuating random flows and compute how this steady state is approached. For 
incompressible flows, the Saffman- Turner expression is an upper bound. 
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1. Introduction 

Turbulent aerosols are of interest in a variety of natural and technological systems. Two 
very important examples are water droplets in turbulent rain clouds [1] and dust grains 
in turbulent accretion disks around growing stars [2j. In both of these systems the 
aerosol is unstable because the suspended particles coUide (leading to aggregation or 
possibly fragmentation of the aerosol particles). The collision processes therefore have 
significant consequences: the formation of rain in one case, and the widely hypothesised 
mechanism for the formation of planets in the other. Collisions always occur due to 
molecular diffusion, but (as pointed out by Smoluchowski [31 H]), macroscopic fluid 
motion can considerably increase the collision rate. 

If the suspended particles are sufficiently heavy (so that their inertia becomes 
relevant), they can move relative to the fluid. In this case the occurrence of 'caustics' 
will typically increase the collision rate by several orders of magnitude O [6]. If the 
aerosol particles are sufficiently light, their 'molecular' diffusion can make a significant 
contribution to the collision rate, which can be estimated by standard kinetic theory. 

In this paper we are concerned with the effect of macroscopic motion of a fluid on 
small particles which have insignificant inertia, so that they follow the flow (advective 
motion). The seminal papers in this area were due to Smoluchowski [3], who flrst 
considered the effect of shear of the fluid flow on collisions, and Saffman & Turner |4j|, 
who gave a formula for the collision rate which has formed the basis for most subsequent 
work on this problem. Their paper was motivated by a problem in meteorology 
argued that small-scale turbulence in convecting clouds can accelerate collisions between 
microscopic water droplets, thus initiating rain formation: in this case the particles are 
brought into contact by hyperbolic or shearing motions of the turbulent flow. 

The formula for the collision rate by Saffman & Turner [1] has been used frequently 
in the past flve decades in cloud physics and in chemical engineering problems. It appears 
to be widely accepted that their expression is an exact relation for the collision rate in a 
dilute suspension. In the following we show that the Saffman- Turner estimate describes 
an initial transient of the problem only. For particles suspended in incompressible 
flows, the collision rate falls below the initial transient (which thus constitutes an 
upper bound). For particles advected in a compressible flows, however, homogeneously 
distributed particles will cluster in a compressible fluid (see for example [3 E]). The 
clustering may increase the collision rate beyond its initial transient. 

The Saffman- Turner approximation treats the flow surrounding a test particle as 
if it were a steady hyperbolic flow, while in reality the flow fluctuates as a function of 
time. In section 13.21 below we give an extension of the Saffman- Turner formula which 
does give the collision rate exactly. Unfortunately the formula contains information 
about the time-dependence of the flow, and it is impossible to evaluate it in the general 
case. 

Because of the importance of understanding collision rates for aerosol particles, it is 
desirable to flnd exactly solvable cases which can be used as a benchmark for numerical 
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studies. The collision rate must depend on a dimensionless parameter describing how 
quickly the fluid velocity fluctuates, the 'Kubo number' [9]. We are able to obtain precise 
asymptotic results on the collision rate in the limit where the Kubo number approaches 
zero. In this case, particle separations undergo a diffusion process. By solving the 
corresponding Fokker-Planck equation we can determine the collision rate exactly. 

The remainder of this paper is organised as follows. In section [2] we introduce 
the equations of motion and the dimensionless parameters of the problem. Section [3] 
discusses the Saffman- Turner theory and our extension of it. Section [3TT] describes the 
expression for the collision rate given in [1] which is the starting point of our discussions. 
(Some new results on the evaluation of the Saffman- Turner expression are described in 
the appendix). In section [3^ we discuss our exact formula for the collision rate, and 
explain why the Saffman- Turner approximation describes an initial transient only. In 
sections HM] we discuss how exact asymptotic results may be found for the limit of small 
Kubo number. A Fokker-Planck equation for the probability density of the separation 
of particles is described in section HI In section [5] this is used to obtain the steady-state 
collision rate and section [6] gives the full time dependence of the collision rate. These 
results are compared to numerical simulations in section [TJ which also contains some 
concluding remarks, discussing scope for further work in this area. 

Finally, we remark that a brief summary of some of the results of this paper has 
already been published [lO] . Here we discuss the problem in greater depth and generality, 
and derive expressions for the time-dependent collision rate which were not discussed in 



2. Equations of motion and dimensionless variables 

We consider spherical particles of radius a in a fluid with velocity fleld u(r,t) which 
has an apparently random motion, usually as a result of turbulence. We assume that 
the suspended particles do not modify the surrounding flow. When the inertia of the 
particles is negligible, they are advected by the flow: 



It is assumed that direct interactions between the particles can be neglected until they 
collide. In other words, the particles follow equation ([1]) until their separation falls below 
2a. 

We model the complex flow of a turbulent fluid by a random velocity field u{r,t). 
We consider flows in both two and three spatial dimensions and for convenience we use a 
Gaussian distributed fleld when we carry through concrete computations. In most cases 
we are concerned with incompressible flow, satisfying V ■ u = 0. Particles floating on 
the surface of a fluid may experience a partly compressible flow [II1IT2], as may particles 
in gases moving with speeds comparable to the speed of sound. For these reasons we 
also consider partially compressible flows. 




(1) 
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It is convenient to construct the random velocity field u{r,t) from scalar stream 
functions or potentials [9J. In two spatial dimensions we write 

u = Af2{V A ijh, + (2) 

where N'2 is a normalisation factor and ip and are independent Gaussian random 
functions. We shall use angle brackets to denote averaging throughout. The fields 
(f){r,t) and tp{r,t) have zero averages, (0) = 0, (ip) = and they both have same 
correlation function, C{R,T): 

(0(r,t)0(r',t')) = C(i?,T) (3) 

where R = \r — r'\ and T = \t — t'\. This two-point correlation function C{R,T) is a 
smooth function decaying (sufficiently rapidly) to zero for large values of R and T. For 
/3 = 0, the flow (I2I) is incompressible. For finite values of /3 it acquires a compressible 
component. In the limit of (3 00 the fiow is purely potential. In some cases the 
physics of a problem dictates that ip and should have different correlation functions; 
many of our results can be generalised in this way. 
In three spatial dimensions we write 

w = A/'3 ( V A A + (4) 

where A = {Ai, A2, A^) and are four independent scalar fields with zero mean and 
with the same correlation function C{R,T); and A/3 is a normalisation factor. 

In the remainder of this paper we choose the normalisation factors to be of the form 

A/". = ^° = (5) 

^-d{d -1 + (3^)C" (0,0) 

where Uq is the standard deviation of the magnitude of the velocity and C" denotes the 
second derivative of the correlation function ([3]) with respect to its first argument. 

Fully-developed turbulent fiows have a power-law spectrum in the inertial range, 
covering a wide band of wavenumbers [131 E]- This feature can be incorporated by 
giving C {R, T) a suitable algebraic behaviour over a range of values of R, as explained 
in [To]. The long-ranged behaviour of the velocity field is not, however, relevant to 
the advective collision mechanism. It therefore suffices to consider a model with a 
short-ranged velocity correlation: for the numerical work reported in this paper we used 
following form of the correlation function 

C(i?,T)=Coexp(-^-^) (6) 

where Cq is a constant. In our numerical simulations we represent the fiow field by its 
Fourier components which are subject to an Ornstein-Uhlenbeck process as suggested 
by Sigurgeirson & Stuart [H]. 

Our problem is characterised by the six parameters listed in table [H From the 
parameters in table (H three independent dimensionless combinations can be formed: 

Ku = mot/?7, riQa'^ , d = 2a/7]. (7) 
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Parameter Symbol 

Particle size a 

Typical velocity fluctuation uq 

Correlation length of the flow r] 

Correlation time of the flow r 

Number density of particles no 

Compressibility /3 

Spatial dimension d 



Table 1. Parameters of the model. 

The flrst parameter characterises the dimensionless speed of the flow and is called Kubo 
number, discussed in [9]. Note that Ku 1 is not possible, because the motion of the 
fluid places an upper limit on the correlation time. Steady, fully developed turbulence 
corresponds to Ku ~ 1. The Kubo number can be small for randomly stirred fluids. It 
is only in the limit Ku — > that we are able to obtain precise and explicit estimates for 
the collision rate: this case is considered in sections HH The second parameter in (JTj) is 
the packing fraction of particles. Throughout we assume that this parameter is small. 
Similarly, the third parameter in ([71) is usually taken to be small. 

We note that Kalda [15] has considered the collision rate for particles in a non- 
smooth velocity fleld: this could be relevant to the case where a ^ rj. 

3. Evaluating the collision rate 
3.1. The Saff man- Turner expression 

Throughout this paper we consider the rate of collision of a given particle with any other 
particle, denoting this by IZ. Some papers consider the total rate of collision per unit 
volume. If the spatial density of particles is no, the total rate of collision per unit volume 
is \nQlZ (the factor of | avoids double- counting). We will not be concerned with what 
happens after particles undergo their flrst collision: in different physical circumstances 
they may coalesce, scatter, or fragment, but in this paper we are concerned only with 
their flrst contact. 

We assume that the particles are spherical (or circular, in two-dimensional 
calculations) and that they all have the same radius, a. We regard the particles as having 
collided when their separation reaches 2a, and we neglect effects due to the interaction of 
the particles and the fluid. In practice the fluid trapped between approaching particles 
may slightly reduce the collision rate , but this effect can be accounted for by replacing 
the radius by an effective radius. 

The problem of calculating the rate of collision therefore reduces to the following 
problem. We consider a given particle, and transform to a frame where the centre of this 
particle is at the origin, and the separation of the centre of another particle is denoted 
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by R. Initially, the reference particle is surrounded by a gas of particles with spatial 
density uq. We assume that these are initially randomly distributed, apart from the 
constraint that none of the particles is in contact with the reference particle. Collisions 
with the test particle occur when other particles come within a radius 2a of the reference 
particle. The rate of collisions is therefore the rate at which particles cross a sphere of 
radius 2a centred at the origin in the relative coordinate system. This is obtained by 
integrating the inward radial velocity over the sphere, and multiplying by the density 
no- This approach gives an expression for the collision rate which we term TZo: 

no = -noj dnvr{2a,n,t)e{-Vri2a,n,t)) . (8) 

where Vr{R,^,t) is the radial velocity at spherical coordinate Q, radius R and time 
t. The function Q{x) is a Heaviside step function, which is used to select regions of 
the surface where the flow is into the sphere of radius 2a. This is the fundamental 
expression for the collision rate given by Saffman & Turner In section [3^2] we discuss 
why this expression is not exact, and give the precise formula. The remainder of this 
section considers how this expression is evaluated; some new results are presented in the 
appendix. 

The evaluation of ([H]) is greatly simplified in the case where the particles are small, 
in the sense that a/r] -C 1. In this case the relative velocity is accurately approximated 
using the velocity gradient of the random velocity field u. This approximation was 
also considered by Saffman & Turner [4j. To lowest order in R, the relative velocity 
R = u{r + R,t) — u{r,t) is approximated as A[r,t)R where A is the rate of strain 
matrix of the flow u, with elements Aij = dui/drj (and i,j = 1, . . . ,d). Two particles 
with radii a <^ t] moving close to each other will thus experience a relative velocity 
AqR. This is illustrated in figure [1] for a flow which is hyperbolic in the vicinity of the 
reference particle. The corresponding relative radial speed Vr is Vr = u^AqR where h 
is the radial unit vector. In particular, at the distance R = 2a where particles collide, 
the radial speed is Vr = 2an^ KqU. Thus when a/rj -C 1, equation (IHD reduces to 



TZo = —2ano I df2 n^Aon 0(— Aon) . 



(9) 



This expression was also obtained in j^. 

The remainder of this subsection is concerned with the evaluation of (Q. Earlier, 
Smoluchowski |3j had considered the special case of a fluid flowing with a uniform shear 
in 3 dimensions: 



(10) 



He obtained the collision rate: 

7^o = ^anoi2af (11) 

which is in agreement with the result of evaluating (Q for this case. For a general 

strain-rate matrix Aq the evaluation of is, however, very difficult. In Appendix A we 
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Figure 1. Schematic picture of two particles of radius a passing each other in a 
hyperboUc flow. The particle at the origin will see particles move past on hyperbolic 
trajectories. Collisions occur whenever particles approach closer than 2a. The collision 
rate is thus determined by the influx of particles into a disc of radius 2a (dashed line) 
around the origin. 



discuss how this expression is evaluated for a general matrix Aq in two dimensions, and 
for a general traceless matrix (representing an incompressible flow) in three dimensions. 

In a turbulent flow, Saffman & Turner [1] argued that the elements of Aq change as 
a function of position and one needs to average over the ensemble of strain matrices Aq 
at different positions in order to estimate the collision rate, so that ([9]) is replaced by: 

7^o = -2ano(^ j dVL fi^Aon 0(-n^Aon)^ . (12) 

At first sight the requirement to average over Aq appears to complicate the problem. 
However, for a rotationally invariant ensemble of random flows, the problem is 
considerably simplified by taking the average. In an incompressible flow, for each 
realisation of u{r, t), the currents into and out of the collision region (the disk or sphere 
of radius 2a) cancel precisely. Therefore ( [T2l) can be written as 

TZq = anQ( / dr2|n^Aon|). (13) 



The same result also holds for cases where the flow is compressible, of the form ([2]) or 
dl]), with /3 7^ 0. This is demonstrated by the following argument (here we discuss the 
two-dimensional case). If we did not include the factor 6(— n^Aon) in f|T3|) . we would 
calculate the sum of the collision rate for the flow u and for the time-reversed flow —u. 
The time reversed flow is generated by reversing the signs of the potentials and ip in 
([2]). The probability density for the Gaussian field —(f) is the same as for (p (and similarly 
for tp). It follows that the collision rate for the time- reversed flow —u is the same as for 
u. The expression f|T3|) is therefore also valid for the cases of compressible flow which 
we consider in this paper. 
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Now using rotational symmetry, one finds the very simple expression 

7^o = anoA,(2a)(|An|) (14) 

where Aa{2a) is the the area of the sphere of radius 2a in d dimensions (explicitly 
A2{r) = Tcr^ and A^i^r) = Airr^/S). For the Gaussian model flow which we consider, 
equation (IT^ gives 

7^o = ^A,{2a)2a^{AnAn) (15) 

V 27T 

with 

iA,,A,,)^-m^i±p^:^. (16) 

3.2. An exact expression for the collision rate 

The Saff man- Turner expression for the collision rate, equation ([8]) or f|T4l) . correctly 
describes the initial collision rate. There are, however, two reasons why the collision 
rate may approach a significantly different value after an initial transient. 

The first reason why ([H]) may fail arises from the fact that the flow field fluctuates 
in time. Recall that we are concerned with the rate at which pairs of particles collide for 
the first time. If the flow is time-dependent, the relative position coordinate R{t) may 
pass through the sphere of radius 2a more than once. This effect may be accounted for 
by writing the collision rate in the form 

n{t) = -n{t) j dnvr{2a,n,t)e{-vr{2a,n,t))x{2a,n,t) . (17) 

As before dQ is the d-dimensional surface element at R = 2a and Vr is the radial 
velocity component (the relative speed) and the Heaviside step function ensures that 
only particles entering the sphere contribute to the collision rate. The factor n{t) is the 
density of particles in the neighbourhood of the test particle at time t. The function x 
is an indicator function: it is equal to unity if the point reaching the surface element 
Q at radius 2a and at time t has not previously passed through the sphere of radius 
2a, otherwise it is zero. The effect of including the function x is illustrated in figure [2J 
The collision rate will reduce below the Saff man- Turner estimate for times where x is 
no longer unity.. 

In the general case (ITTl) is difficult to evaluate since the indicator function x depends 
on the history of the flow. If the flow is rapidly fluctuating (that is, if the Kubo number 
of the flow is small), the relative separation of two particles undergoes a diffusion process 
which makes it possible to exactly evaluate the collision rate for a given ensemble of 
random flows. 

A second effect may modify the collision rate from the Saffman- Turner estimate 
is that if the flow field is compressible, particles may cluster together, and the particle 
density in the vicinity of a test particle may be higher than the average density. This 
effect is expected to cause the collision rate to increase, after an initial transient during 
which the clustering becomes established. Again, this effect is very hard to quantify in 
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Figure 2. Initial positions of particles which have collided with a reference particle at 
the origin after a certain time (shaded regions) . The collision region (disk of radius 2a) 
is bounded by a dashed circle. At short times this region extends along the direction 
of the stable eigenvector of the initial flow. At intermediate times, we see the effect of 
this eigenvector having rotated. At long times, the set of colliding initial conditions is 
stretched and folded. The shaded region then covers most of the surface of the circle, 
including the region where the radial velocity is inward. The figure shows results of a 
simulated flow with small Kubo number, of the type discussed in section 21 at times 
t = 20r, t = lOOr and t = lOOOr. 



the general case, but we can obtain precise asymptotic results in the limit of small Kubo 
number, where a diffusion approximation becomes applicable. Because we consider 
particles which are advected by the flow, there is no clustering effect when the flow is 
incompressible. 

We note that for incompressible flow the only correction to the Saffman- Turner 
estimate is the occurrence of the factor x in dill)- This implies that for advected particles 
in an incompressible flow, the Saffman- Turner estimate of the collision rate is an upper 
bound. 

4. Diffusion approximation for small Kubo number 

In the limit of a rapidly changing flow, that is when Ku ^ 1, particle separations 
R = r—r' undergo a diffusion process (see [8] for a review). By solving the corresponding 
Fokker-Planck equation with the appropriate boundary conditions, we can determine 
the collision rate exactly in this limit. 

4-1. Fokker-Planck equation 

Consider two particles, one at r and the other in its vicinity, at r'. The equation of 
motion ([T]) implies that their separation R = r' — r obeys 





Integrating (fTSjl over a short time interval 6t we obtain 




(19) 



10 



The moments of the components SRi of SR are: 

/t+5t 
dt' {u,ir + R,t')-u^ir,t')) = 

/oo 
dt{{ui{R,t) -Ui{0,t)){uj{R,Q) -Uj{0,0))) , (20) 
-oo 

assuming that 6t ^ r. We thus obtain the following Fokker-Planck equation for the 
density p{R,t) of particle separations fl\ [8l [TO]: 

|^ = V^D(H)Vp. (21) 

Here 3{R) is a diffusion matrix with diffusion coefficients 

D,,iR) = ^^{6R,6R,) 

/oo 
dt [{u^iO, t)ujiO, 0)) - {uiiR, t)ujiO, 0))] . (22) 
■oo 

In terms of the correlation function C {R, T) introduced in section [2] we have 

D,,{R) = -Uj JZ, dt [6,, ((rf - 1 + (3')C"iO, t) 

^±mlc'{R, t) - C"{R, t)) (23) 

+{l-l3^)^{-lC'{R,t) + C"{R,t)y . 

In order to determine the collision rate it is convenient to write the Fokker-Planck 
equation in the form of a continuity equation 

^ + V^j = (24) 

where the components of the probability current can be identified as j = —3{R)'V p. 
The collision rate 7?.(t) between particles of radius a is given by the rate at which the 
particle separation R decreases below 2a ^Oj. This rate equals the radial probability 
current of particle separations evaluated at i? = 2a, i.e. 

7^(^) = - J dn fi^j{R = 2a, t) (25) 

where n is the radial unit vector. We evaluate fl2B]) in two steps. First (in section 
we determine the steady-state collision rate as t — oo. Second, we determine the full 
time-dependence of the collision rate (section E]), describing how the initial transient 
discussed in section |3] approaches the steady state. 

4.2. Transformation to spherical coordinates 

Because of the angular symmetry of the statistics of the fluid velocity, it is convenient 
to transform (12T!) to spherical coordinates. The probability density p depends on the 
separation R only and obeys 

R'-'^{R,t) = -^{R'~'UR,t)) . (26) 



11 



with 



Here the radial probabihty current jr{R,t) is given by 

jr{R,t) = -(3^U2f{R)p{R,t)+U2g{R)^p{R,t) , (27) 

f{R) = dt[i^ (^c'{R,t)-C"{R,t)^ +C"\R,t) , 
g{R) = dt {d-l+ f3^)C"{0, t) - ——C'{R, t) - (3^C"{R, t) . (28) 

J-oo L -I 

5. Steady-state solution 

Consider now the steady-state solution of ( l2T|) . obeying dp/dt = 0. The steady-state 
collision rate, denoted by TZoo, is determined by (125|) in terms of the steady-state current. 
The total current —Ad{R)jr{R) entering a sphere of radius R and area Ad{R) must be a 
constant. This gives jr{R) = —T^oo/A.d{R) where T^oo is a constant, equal to the collision 
rate which we wish to determine. 

The equation of motion ([T]) is only physically meaningful when a -C 77, so that the 
fluid velocity is approximately the same throughout all of the region occupied by the 
particle. We will, however, solve the diffusion equation for general a. One reason for 
treating the general case is that it is hard to verify the analytical expressions for the 
limiting case 0/77 -C 1 in numerical investigations. 

To determine this current, it is necessary to consider the appropriate boundary 
conditions. First, if the particle distribution initially is uniform with density no, we 
expect it to remain so at large separations. Thus we use the boundary condition 

p{R-^ 00, t) = no. (29) 

The boundary condition ( !29l) is implemented as follows. For a specified current density 
jr{R,t) we can solve ( 1271) by finding an integrating factor h: 

p(fi.t)^Mfi)(A-/;d/?j^) (30) 

with integration constant A and 



00 



h{R) = exp ( / di?'4^ ) • (31) 



R 



9{R') 

The boundary condition fl29l) determines the integration constant in fl30|) to be Aq. 
Second, when particles comes closer than the distance of 2a they collide and must be 
removed. This is taken care of by the following boundary condition at 2a: 

p(2a,t) = 0. (32) 

Inserting this boundary condition into (1301) yields 
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Using jr{R) = —7loo/Ad{R') and solving for the constant T^oo we find the following 
expression for the collision rate: 

= -rm^' [L '''5(5)'="= r L \ 

where T{z) is the Gamma function. Equation (1341) is the main result of this section. We 
continue by discussing a number of limiting cases. 



5.1. Collision rate for incompressible flows 



For an incompressible fiow (/5 = 0), the integrating factor is just unity (incompressibility 
'V'^u = implies V"^© = which in turn implies that the term proportional to p in 

simplifies to 



(127|1 vanishes). In this case, equation 



uquI 



T{d/2) dC"{0,0) 



2a 



dRU dt(^R'^-^C"{0,t) - R'^~^C'{R,t) 



. (35) 



5.2. Collision rate for f3 



When (3=1, the strength of the solenoidal and potential parts of the field are equal. This 
case may be relevant to the dynamics of particles floating on the surface of a turbulent 
fluid The collision rate can be evaluated exactly in this case, by rewriting ( IMl) as 



7^n 



,/32 Ing(oo) 



T{d/2)d{d-l + (3^)C" (0,0) 

dRR^-'^g{R)^'-^exp(/3\/3^-l) / dR 



2a 



J ZdtC'"{R',t) 
9iR') 



-1 
(36) 



Setting = 1 we flnd 



7^. 



UquI 



T{d/2)d^C"{0,0) 

J°l dt {dC'iO, t) - lim«,_„o '^C'iR!, t) - C"(oo, t)) 



X 



jZ<^RR'- 



(37) 



for d < 2. When d > 2 we obtain 



7^n 



(rf-2)(2a 



271^^/2 



X 



r(rf/2)rf2C"(o,o) 
d 



dt rfC"(0,t)- lim ——C'{R\t)-C"{oo,t) 



(38) 



R'~*oo R' 

What makes it possible to flnd exact solutions in this case is the fact that when (3 = 1, 
the functions f{R) and g{R) in fl28l) are related as g'{R) = f{R)- Note that this relation 
is always true in one spatial dimension, independent of the value of (3. Note also that 
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when g'{R) = f{R), equation (ITTI) can be solved for p by integration (starting at e.g. 
2a, with p{2a,t) = 0) 

p^R) = ' / dHjiB!) = / dR'- 



9iR)J2a ' 9{R)]2a Aa{R!) 

^_r(rf/2)^ rdfi"^ 

27r'^/2 g{R) J,^ R"i- 

Thus, a non- vanishing steady-state colhsion rate is obtained only if g'{R) 7^ f{R) or 
when d > 2. 



• (39) 



5.3. Collision rate for small particles 

The collision rate must vanish as the particle size approaches zero. This means that the 
integral in the denominator in must diverge for small a. Thus if a -C is small 
enough, the major contribution in the integral in comes from small values of R and 
the relative error will be small if we replace the integrand by its small R expansion. It 
is convenient to change variables according to R = R/rj. We expand 

C{R) = C(0, T) + ^dlC{0, T)R' + ^4^(0, T)R' + ■■■ (40) 
and obtain approximately 

- (3^Mlf{R) = -{d + l + {l~ d)T)Vr]R +■■■ 
Mh{R) =-V{^Rf + ... 

h{R) = R-id+^+i^-'i)^) + . . . . (41) 

Here the parameter T> is defined as 



^ 1 d^Dii 
T) = - 



2 di?2 



i?i^O,Rj^i=0 

d + 2 ulr]-^ 



dt9|C(0,t). (42) 



3\d{l + {d - 1)T) dlC{0,0) 

It corresponds to the diffusion constant Di/(m7)^ in equation (39) of [9] and in the 
incompressible case it corresponds to the diffusion constant T> defined in equation (17) 
of [To] : here F is an alternative parametrisation of the compressibility, defined as follows 

^ ^ D22{Rl^0,R^^l = 0) ^ d+l+[3' 

The parameter F was also used in [H] (and earlier work cited therein) . It can take values 
between F„ji„ = 1/3 and F^^^ = {d+ l)/{d — 1). Here T^^^^ corresponds to a completely 
incompressible flow with (3 = and F,,,;,, corresponds to a purely potential (/5 00) 
compressible flow. In particular, when F = 1 the field strengths of the compressible and 
incompressible components of the flow are equal, /3 = 1. 

Using equations fBTl) . we get the collision rate for small particles 

271^^/2 . /2a\^'^"^^^~^ 
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This expression was previously derived in pTO]. 

For an incompressible flow, with F = {d + 1)/ (d—l), the particle density is uniform 
and the collision rate is proportional to the packing fraction nQa'^ as expected. By 
contrast, if the flow is compressible, the particles will cluster on a fractal set This 
is expected to enhance the collision rate because the particle density is large within the 
clusters. The clustering effect can be characterised by the correlation dimension of the 
particles, D2. It is deflned by the scaling law P(e) oc e^^, where P(e) is the probability 
that the particle separation is smaller than e p2]. We have 

P(e) oc dR R"-^ p{R) oc e^"^"^)^"^ (45) 
Jo 

and thus D2 = {d — 1)T — 1. This expression was derived in [7j. 

Equation shows that in compressible flows, the collision rate depends upon 
the correlation dimension D2 = (rf — 1)F — 1 < d. When d <^ 1 the collision rate 
in a compressible flow is therefore much larger than the corresponding rate in an 
incompressible flow. 

Note flnally that the steady-state collision rate (jS]) tends to zero when F — >■ 
l/{d — 1), i.e. when /? d{d — l)/{A — d) = /3c- For values of (3 larger than (3c, no 
steady-state current satisfying our boundary conditions exists. In one spatial dimension 
the parameter f3c vanishes: no non-trivial steady state exists because all particle 
trajectories eventually coalesce (this effect was termed 'path-coalescence' Wilkinson & 
Mehlig [IB]). 



5.4- Collision rate for the Gaussian correlation function 

When the correlation function C{R, T) is of the Gaussian form (Q, we obtain from (E 

47r''/2no?7°' 



d{d 
^2 



X exp ( (3 



-P')T{d/2) 
dRR'-'^e^"/^ 

dR'R' 



2 + d-R'^ " 



R 



(rf - 1 + /?2) (eR'V2 _ 1) + /?2^'2 y J J (46) 

where as before R = rjR. The major contribution to the integral comes from small 
values of R (except for when d = 2 and (3 = 1 or d = 1, when the integral diverges) and 
for small values of a, the integrand can be expanded in powers of R. Expanding the 
exponential function in the integrand yields 

2((rf- 1)F- l)7r'^/2. 



7^. 



T{dl2) 

where the function h was deflned in ( I3ip . The diffusion constant V is given by 
With ([6D we flnd 

rf-l + 3/?2 uIt 



(47) 



V 



d{d-l + (3^) r/2 ■ 



(48) 
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For low dimensions, 

h{a) ^ a^'i-i)r~d-i ^^g^ 

is a fairly good approximation (a maximum of 3 percent error when d = 2). Substituting 
(149|1 into (147|) gives once more 



6. Time-dependent solution 

So far we have derived the probability densities and collision rates for the steady state. 
We now wish to derive expressions of these quantities as functions of time. To this end 
we need to solve the time dependent Fokker-Planck equation (!26l) with the boundary 
conditions 

p{2a,t) = , and p{oo,t) = Uq . (50) 
As before we assume a uniform initial scatter of particles 

p{R,0) = e{R-2a)no. (51) 
The collision rate is given by equation ( l25i) 



nt) = fW^>(^ = t) , (52) 

where we have used that the area of a d-dimensional sphere of radius r is Ad{r) = 
27r^/V^-7r(d/2), and the gamma function is denoted by T(z). 

In order to render the boundary conditions homogeneous, we split p{R, t) = 
pi{R,t) + p2{R) and impose the conditions pi(2a,t) = pi(oo,t) = 0, and that the 
differential equation homogeneous in p is also homogeneous in pi. Thus p2{R) is uniquely 
given by the steady-state solutions found in the previous section. 

Now consider the remaining equation for pi{R,t) which is identical to ( !26l) with 
homogeneous boundary conditions and initial condition 

Pi{R,0) = e{R- 2a) {no- p{R, oo)) . (53) 

Separation of variables pi{R,t) = J-'{R)Q{t) in gives (using the radial current ([2ZD) 

^^^g_ _ 1 d 
g~dt ~ R^^dR 

Since the left-hand side depends on t only, and the right-hand side depends on R only, 
both sides must be equal to a constant, B say. Choosing B = —Vp^^ where /i is a 
positive dimensionless constant, we obtain 

Q{t)=goe-^"'^\ (55) 

Note that the solution corresponding to jU = has already been taken care of in p2- 
Therefore we can restrict ourselves to considering p > in the following. 



f{R):F + g{R)j^:F 



(54) 
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To solve the radial part of fl54l) . we consider the limit of i? <^ 77, where f{R) and 
g{R) are simple power laws (this follows from fHT]) ). Using the dimensionless variable 
R = R/t] gives the following approximate equation for small values of R 
_ cfT _BT 

R^~ + (2c/ + 2 + (1 - d)T) R-= + (c/ + 1 + (1 - d)T)dr = -/i2^.(56) 
oR^ oR 

This is an Euler equation. Its solution is obtained by the variable substitution R = ae", 

where the factor a = 2a/rj is included in for later convenience and < u < 00. We find 

I sinh ( fJ?\n{R) + 0+] , if < /i^ < /ig 

where 

1 . . ,s ^ rf(l - rf) + (4 - , , 

^, = -(i + (i-d)r)^ 2(^-1 laW (58) 

and and 0± are integration constants. 

For all values of fiQ in the allowed range — ^ < /io < |(4 — c?), the boundary 
condition jF(oo) = is automatically fulfilled. The boundary condition jF(a) = gives 
(j)± = — A/±(/io — /i^) ln(a) + i/^nvr, and thus 

- ^ I sin ( - /io ln(^/ a) ) , if /^^ > /io 

I D+ sinh f A//io - /i^ ln(i?/a) j , if < /i^ < /ip 

The boundary conditions do not constrain the eigenvalue fi and we need to consider a 
continuous superposition of eigenf unctions 

Pi{R,t) = R-'^-'"" d/iC(/^)sin (^^/i2 -/i2ln(i?/a)^ e"'^'^* 

+ d/i C(/i) sinh (^^fi4^ \n{R/d)^ e"'^'^* . (60) 

Here ({fi) and C(Ai) are functions to be determined by the initial conditions. It turns 
out that it is sufficient to consider ({fi) for our initial condition and take ({^) = from 
now on. 

At time t = we make the change of variables u = ln(i?/a) as before. We set 
s = \/ IJ''^ — iJ-o and find 

p^{u, 0) = — (ae")-'^-^" r ds ^C(s) sin (su) . (61) 
du Jq ds 

Comparing this to the required initial distribution 

Pi(m,0) = ^e(iz)(no-p(M,oo)) (62) 



gives 



[ ds^C(s)sin(sM) = (ae")'^+^'°e(M)(no-p(M,oo)). (63) 
Jo 
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The left hand side is of the form of a Fourier sine transform with inverse 



-((s) = -a'^+'^« / dM(no -p(M,oo))e^''+^°^"sin(su). (64) 

ds TT Jq 

We insert this expression into the ansatz for pi . Upon changing order of the integrations 
(that is we perform the s-integral first), we find 

p^fu t) = ^ ^~iM?,m-(d+„o)u~uy(4m) 

X du (no - p(n, oo))e('^+^«)"-"'/(^^*) sinh (^^^ . (65) 

To calculate the collision rate we need to know the current jV at radius 2a. To this end 
we just require the derivative of p evaluated aX R = 2a (since p itself is constrained to 
vanish there). We find 

(2a, t) = \ / du (no - p{u, oo)) Me('^+^°)"-"'/(^^*) (66) 



Equation (1521) now allows us to calculate the time-dependent collision rate. We split 
ir [see equation fl271) ] into two parts depending on pi and p2 respectively. The second 
part gives the steady-state collision rate T^oo, see equation flM|) . We find (expanding in 
powers of R) 

^{d-i)/2^2aYe-^'^' 



n{t) = 7^c 



r(rf/2) y/V^ 

) 

du (no - p(u, oo)) Me('^+^«)"-"'/(^^*) . (67) 



This is our main result for the time-dependent collision rate, valid for correlations with 
non- vanishing fourth order derivative at i? = and for small particles. 

An approximation to this result can be obtained by approximating the steady-state 
probability density 

m oo) ^ / d^'^^^/ 1 '^R'^^l^^ (68) 

by expanding in powers of R. Since the w-integral extends to infinity and since the 
Gaussian contribution in the integrand becomes small for large values of Vt, we must 
have that p{u, oo) goes to no as u goes to oo, which is the case for the original p{u, oo). 
To accomplish this, we match the small-^ expression of p{u, oo) at ^ = 1 to the large-^ 
expression. By using (jUj) and (jSj) and matching at ^ = 1 we find 

p{R, oo) = nov {Q{R - 1) + R^'^'^^'il - e(i? - 1)) - a-2Mo/^-d) ^ (69) 

where the last term is not cut off at ^ = 1 as it vanishes sufficiently quickly. 

The time-dependent collision rate can now be evaluated from the sum of three 
integrals of the type 
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X erf((6„ + rf + /io)vI't) - erf (6„ + (i + /io)Vl't 



where a„ = {1, -a-'^-^w, a-««-2Mo}^ 5^ = {0,-d - 2/io, -rf} and c, 
{— ln(a), — ln(a), 00}. Putting everything together we find 



(70) 



7^(t) 



2. ^ . ^d. 



r(rf/2) 
1 



e-'^o^* + (rf + /io)e 



1 + erf 

d{d+2^io)Vt 



-fioWt 



In a 



2VVt 



X < erf 



{d + fio)VVt 



— erf 



In a 



2VVt 



(71) 



vahd for a <^ 1. For times less than [d{d + 2fiQ)V] ^ the above expression simphfies to 



2/io 



(72) 



For intermediate times, [d{d + 2fiQ)V] ^ ^ t ^ ln(a)/(2/ioI^), we find approximately 



2flQ 



X < erf 



(rf + /io)Vr't 



erf 



{d + fioWVt + 



In a 



2Vr't 



while for large times, t ^ \n{a)/(2fioD), where 



7^(t) ~ -7^oo 



1 + erf 



In a 



2VVt 



(73) 



(74) 



Consider finally evaluating fl7T|) in the incompressible limit, we have /xq = —d/2 and the 
terms involving the integral cutoffs ln(a) cancel (/i = —I2 in flTUl) ). Thus in the case 
/3 = the collision rate simplifies to 



7^(^) = d 



7i'^/\2ay 

'WJW 



■noVl 1 + 



d 



(75) 



dVirVi 

Equations ( !7T|) and ( !75l) are compared to results of numerical simulations in the following 
section. 



7. Numerical illustration and concluding remarks 

We performed simulations of the collision rate as a function of time at small Kubo 
numbers, for both incompressible and compressible fiows. These illustrate (figure [3l) the 
very complex behaviour of the model. For example, in the case of compressible fiows 
the collision rate at first decreases below the value determined by the Saffman- Turner 
approximation due to the effect illustrated in figure El before rising again as particle 
clustering becomes apparent. 
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Figure 3. Collision rate for particles advected in a two-dimensional flow of the 
form ^ with ^ and Cq = u^rf' /2. Particles are initially randomly distributed, 
initially overlapping particle pairs are removed, as are particles which have collided. 
Parameters: rip = 1000, Uq — \, r] — 0.1 and a = 0.001 for all graphs and for a /3 = 
and r = 0.004, b /? = and r = 0.0008, c (3 = 1/^13 and r = 0.004, d /3 = l/VTS and 
r = 0.0008, e (3 = l/%/5 and r = 0.0008, i (3 ^ .Jzjl and r = 0.0008. The collision 
rate is approximated by the cumulative sum of all collisions up to i, divided by t. The 
light blue areas are the intervals R±2u / ^fN , where a is the standard deviation of the 
rate, and N is the number of realisations of the flow. The latter decreases over time 
(it is typically N ^ IQ^ for the first two decades in t, and ~ 10 for the last decade). 
The Saffman- Turner estimates (|15ll6p are shown as red dashed lines. Our own result 
([7T|) is shown as red solid lines. To correspond to the simulated collision rate, the 
plotted theoretical collision rates has been integrated to and then divided by t. Due to 
the assumptions in the Fokker-Planck theory, the long time rate dTlT) is not valid for 
t < T. The plotted long time theory is a combination of the short time rate ([76l) and 
the long time rate (|7ip . matched at the time at which the long time rate drops below 
the short time rate (which is an upper bound). 
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7.1. Simulations 

Figure [3] summarises our results for the collision rate of particles advected in flows with 
small Kubo numbers. Shown are numerical simulations for a two-dimensional random 
flow of the form ([2]) with correlation function ([6]). 

Panels a and b show the collision rate in an incompressible flow {(3 = 0). As 
expected (see section [3^ . the collision rates drops below the initial transient given by 
the Saffman- Turner approximation, ( fTSfTG]) with /3 = 0. For the particular choice ([6]) 
with Co = u^rf equations fll5lll6p become in two spatial dimensions 



Also shown is our own theory valid for small Kubo numbers (equation (17T1) in section [6]). 
The agreement between the theory and the simulations is good in all cases, but slightly 
better for the smaller Kubo number (panel b). As discussed in section |3T2| the initial 
transient constitutes an upper bound to the collision rate. 

Panels c to f in figure [3] show the collision rate in compressible flow. Now, the 
Saffman & Turner theory is no longer an upper bound (see for example figure [3]F), 
because initially homogeneously distributed particles in an incompressible flow cluster 
together. The corresponding density fluctuations increase the collision rate, as our exact 
result shows (equation ( |7T1) in section [6]). Again we observe good agreement between 
the simulations and our analytical result. 

7.2. Scope for further investigations 

In this paper we have concentrated upon the solvable case of advective collisions in flows 
with small Kubo number, which provides considerable physical insight. We conclude 
by commenting on the relation between these results and collisions in a turbulent 
flow field which satisfies the Navier-Stokes equation. The standard approach, based 
upon the Saffman- Turner formula, predicts a collision rate TZ ~ noa'^T"^, where r is 
the Kolmogorov timescale [13] of the turbulent flow. The calculation based upon the 
diffusion equation gives a collision rate TZ ~ Kunga'^r^^. Given that Ku = 0(1) for a 
turbulent flow, we see that the Saffman- Turner and diffusive expressions are of the same 
order. 

These observations are consistent with the hypthesis that the collision rate for small 
particles in a turbulent flow is 



where 8 is the rate of dissipation per unit mass, v is the kinematic viscosity, and K^i is a 
universal constant (depending only on the dimension). It would be a valuable addition 
to the literature on aerosols and suspended particles to determine the value of from 
numerical simulations using a Navier-Stokes flow. 




(76) 
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Appendix A. Evaluation of collision rate for constant shear 

In this appendix we show how to evaluate the expression ([9]) for a general matrix A 
with elements (we drop the subscript 0). It is convenient to decompose A into a 
symmetric part A() and an antisymmetric part A[] 

Au \ 

421 

2 



^12+^21 4 
^ ^22 



M2-A2, 



+ I -^12^ 



(A.l) 



The collision rate is independent of the antisymmetric part Aq (because rotations do 
not contribute to the collision rate). The symmetric part A() can be diagonalised by 
an orthogonal transformation, A() = OSO"^, where S is diagonal with the eigenvalues 
o"i < o"2 cTn of A(). In evaluating ([9]) we may write n An = n OSO h = n Sn , 

where n' = O^n. Since n' is just a rotation of h which is integrated over all directions, 
we can replace ([9]) by 

7^o = -2ano j d^] n^Sn 0(-n§n) . (A.2) 
Appendix A.l. Two spatial dimensions 

We now show how to perform the integral in (1A.2P in two spatial dimensions. Note that 
if the flow determined by A is not area preserving, the particle density will change as a 
function of time. In this case uq in ([9]) must be replaced by 

n(t) = no^^ = noexp[-trAi(:] . (A. 3) 

At short times we approximate n{t) ~ uq and find 

^ -{2aynQ / d^ (cxi cos^ 9 + a2 sin^ 9) e(-ai cos^ 9 - a2 sin^ 9) 
Jo 

—7ia^ if cr+ < — (T_ 

2(2a)^no < a/ct^ — o-+ — (7+ arccos(o-+/(j_) if — a_ < < a_ (A. 4) 
if cr+ > (T_ 

where a+ = (ai +(T2)/2 and cr_ = (0-2 — cri)/2 > 0. This is the final result in two spatial 
dimensions, expressed in terms of the eigenvalues ai < (T2 of the symmetric part A() of 
the strain matrix. 



Appendix A.2. Three spatial dimensions 

In three spatial dimensions we only consider incompressible shear fiows, that is a general, 
three-dimensional traceless matrix Aq. It has eigenvalues o"3 > cr2 > ai, obeying the 
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relations 0-3 > 0, —0-3/2 < (T2 < (73 and ai = — (T2 — (73. We use a spherical coordinate 
system, h — (sin 9 cos (f, sin ^ sin (/?, cos 9) and find 

n^Sfi = (—(72 — (73) sin^ ^ cos^ </? + (72 sin^ ^ sin^ </? + (73 cos^ ^ 

^ (73 I ^ cos^ ^-l-(^^ + 0(l - cos^ ^) cos(2^) I . (A.5) 

Substitute t — cos 9 and 2(p ^ (p^ — (p we obtain 

7^o = ((72, (T3) = - / dv? / d^sin(^)n^§ne(-n^Sn) 
Jo Jo 



- 2 / d(^ / 



dtn'Sne(-n'Sn), (A.6) 

Jo Jo 
where the transformed integrand 

n^Sn = ^3 { - ^ - + (1 - t^) cos ^} (A.7) 

is smaller than for t < to, where to is given by 

9 1 + X cos 09 , . 

^0 = 77- A.8 

3 + X cos LP 

and a; = (2(T2/(T3 + 1). Because the values (T2 can take are limited by (T3, x must be in 
the range < a; < 3 and to in the range —00 < to < \- Performing the t-integral from 
to to gives 

. / -t3 + to + a;(to-|)cos<^ ifto>0 

^0(0-2, (73)= / a(/7(73< V 

Jo [0 if to < 

2 r° , (1 + xcos<^)^/^ 

= 3^3/ ^ (A.9) 

-J J -m (<J + X cos 99) 

where (/Pq = arccos(— l/x), if (72 > and (p^ = tt, if a2 < 0. 

Expanding the integrand around in x around the point Xq (to be determined below) 
we find 

00 n 

72 = 2(73 ^ ^ (^) cos" (^(-2)-"(2n -2k- 5)\\{2k - 1)!! 

n=0 fe=0 

X (1 + Xo COS (3 + Xo COS (^^^ _ (^ -^Q^ 

77/. 

If we choose xq = 0, which corresponds to (T2 = —(73/2 < 0, we can perform the 
integration after changing the order of summation and integration. The collision rate 
becomes 

no{a2<0,a3) (A.ll) 

_47r ^^ {An-2k- 5)!!(2fc - l)!!(2n - 1)!! /(72 iV" 

'V^'^'hh 2"3'=(2n-fc)!fc!n! 2 J " 

Only even powers contribute to this sum: we have not used that o"2 > and could thus 
as well have expanded the starting equations using ai and (T3, with the only difference 
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that (72 would be replaced by —cr^ — cti in (lA.lip . We thus find that 7^o(<^2)0'3) is 
symmetric around ai = (72 = — (73/2, i.e. 7lo{a2, cts) = 7lo{—a2 — Cs, cts). 

To obtain an expression valid for a2 > 0, we reorder the eigenvalues of the matrix 
A() as (T3 — » cTi and cxi — (T3 = — CT2 — giving the eigenvalue ranges cti < and 
o'l < 0-2 < —ai/2. The integrand analogous to ( 1A.5I) becomes 



(A.12) 



where cxi has opposite sign as 173 before. This integrand is smaller than for t > t^, 
where to is given by 

ti = l±y^ (A.13) 

° 3 + ycos(/? ^ ^ 

where y = (2(T2/cri + 1) lies in the interval < y < 3 and to lies in —00 < tl < 1/2. 
Performing the t integral from to to 1 gives 



TZo{a2,(T3] 



2n 



dipa 



3?/ cos ip + tQ — to + y{-^—toj cos ip if to > 



\3/2 



3?/ COS ip 

- (1 + w COS 09)" 

dv9 Yn 

■ipo (3 + ?/ cos ip) 



if tn < 



(A.14) 



where we have proceeded as in deriving flA.12p with ipQ = arccos(— l/y), if 0-2 < and 
(fo = n, ii a2 > 0. We obtain the final result 

2n 



no{a2,a3) = no{2af< 



ra=0 

47r(o-2+o'3) ' 
v/3 



(A.15) 



n=0 



where 



2n 

E 

k=0 



(4n-2fc-5)!!(2fc-l)!!(2n-l)!! 
2"3'=(2n - k)\k\n\ 



(A.16) 



This is the final result for incompressible flows in three spatial dimensions. 
For the particular case ( fTOll considered by Smoluchowski [3] we have 

/ a/2 \ 



A 







(A.17) 




y a/2 J 

with eigenvalues —a/2,0, a/2. Substituting (72 = and = a/2 into (1A.15I) we 
obtain ffTTl) . When a2 = this expression agrees with the classical result ffTTl) due to 
Smoluchowski [3J. 



